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Linear Paul traps have been used recently to simulate the transverse field Ising model with long- 
range spin-spin couplings. We study the intrinsic effects of phonon creation (from the initial phonon 
ground state) on the spin-state probability and spin entanglement for such quantum spin simulators. 
While it has often been assumed that phonon effects are benign because they play no role in the 
pure Ising model, they can play a significant role when a transverse field is added to the model. 
We use a many-body factorization of the quantum time-evolution operator of the system, adiabatic 
perturbation theory and exact numerical integration of the Schrodinger equation in a truncated 
spin-phonon Hilbert space followed by a tracing out of the phonon degrees of freedom to study this 
problem. We find that moderate phonon creation often makes the probabilities of different spin states 
behave differently from the static spin Hamiltonian. In circumstances in which phonon creation is 
minor, the spin dynamics state probabilities converge to the static spin Hamiltonian prediction at 
the cost of reducing the spin entanglement. We show how phonon creation can severely impede the 
observation of kink transitions in frustrated spin systems when the number of ions increases. Many 
of our results also have implications for quantum simulation in a Penning trap. 
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I. INTRODUCTION 

Complex states of matter like spin liquids are suspected 
to exist in quantum spin models with frustration due to 
geometry or due to the nature of the spin-spin interac- 
tion 11-3 . Spin liquids are complicated quantum many- 
body states that exhibit significant entanglement of their 
wave functions without symmetry breaking, and could 
also exhibit emergent quantum phenomena within their 
low-energy excitation spectra. Classical computation, 
such as exact diagonization and quantum Monte Carlo 
simulation, or conventional theories based on local order 
parameters fail to describe these systems without bias. 
For example, exact diagonalization studies are limited to 
small size lattices and hence usually have strong finite- 
size effects, while quantum Monte Carlo simulations can 
suffer from the sign problem or have a large computa- 
tional expense to describe long-range interactions and 
hence cannot reach the low temperatures needed to see 
the predicted exotic phases. 

Feynman proposed that one could use controlled 
quantum-mechanical systems with few quantum gates to 
simulate many-body problems [HIS] as an useful quantum 
computation before achieving universal quantum com- 
putation. In recent years, there has been significant 
success in trying to achieve this goal by quantum sim- 
ulation of desired spin models through analogous cold 
atom systems [BHH]- We focus here on one platform 
for performing analog quantum computation, the sim- 
ulation of interacting quantum spins via manipulation 
of hyperfine states of ions in a linear Paul trap [51115] 
although many ideas presented here can be generalized 
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to adiabatic quantum state computation in the two di- 
mensional Penning trap as well |5]. In the Paul trap 
systems, clock states of the ions (states with no net z- 
component of angular momentum) are the pseudospin 
states, which can be manipulated independently by a 
pseudospin-dependent force driven by laser beams. The 
lasers couple the pseudospin states to the lattice vibra- 
tions of the trapped ions, which leads to effective spin- 
spin interactions when the phonon degrees of freedom are 
adiabatically eliminated [TUl HU [H] based on the idea of 
geometric phase gate [12] or M0lmer-S0rensen gate |17) . 
Theoretically, the analog ion-trap simulators can be de- 
scribed as nonequilibrium driven quantum systems with 
both spin and phonon degrees of freedom. Sufficiently 
small systems can be treated numerically in an exact 
fashion by truncating the phonon basis and taking into 
account all possible quantum states in the solution of the 
time-dependent Schrodinger equation. Experimentally, 
ion traps have been used to simulate the transverse-field 
Ising model with a small number of ions [HI HH dl] based 
on simulated quantum annealing [TH] (see Ref. [121 for a 
review). It has been known experimentally that moder- 
ate phonon creation is commonplace (on the order of one 
phonon per mode) |1T] . even when the system is cooled 
to essentially the phonon ground state prior to the start 
of the simulation. In addition, the role phonons play are 
intrinsic and essential for the mediated spin-spin inter- 
action in trapped ion systems especially in the presence 
of noncommuting magnetic field Hamiltonian in addition 
to the spin Hamiltonian of interest. Therefore, an under- 
standing of the role phonons play in the spin simulator 
is crucial to understanding its accuracy. 

The organization of this paper is as follows. In Sec. 
II, we describe the microscopic Hamiltonian for the ion- 
trap-based simulators and then show how one can factor- 
ize the time-evolution operator into a pure phonon term. 
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a coupled spin-phonon term, a pure spin-spin interaction 
term, and a complicated term that primarily determines 
the degree of entanglement of the spins. Next we use adi- 
abatic perturbation theory to determine how adiabatic 
state evolution can be used to reach a complicated, po- 
tentially spin-liquid-like ground state, and detail under 
what circumstances the evolution is not adiabatic (dia- 
batic). In Sec. Ill, we show numerical comparison stud- 
ies in various relevant circumstances based on a direct 
integration of the time-dependent Schrodinger equation, 
including both spin and phonon degrees of freedom (the 
latter in a truncated basis). In Sec. IV, we conclude with 
discussions and possible experimental limitations and im- 
provements. 



II. THEORY 
A. Microscopic Hamiltonian 

When N ions are placed in a linear Paul trap [^ fTUll^ 
with harmonic trapping potentials, they form a nonuni- 
form (Wigner) lattice, with increasing interparticle spac- 
ing as one moves from the center to the edge of the 
chain. The ions vibrate in all three spatial dimensions 
about these equilibrium positions [21] with 3A^ normal 
modes. Two hyperfine clock states (relatively insensi- 
tive to external magnetic field fluctuations because the 
z-component of total angular momentum is zero) in each 
ion will be the pseudospins (and are split by an energy 
difference huJo). Hence, the bare Hamiltonian Hq includ- 
ing the pseudospin and motional degrees of freedom for 
the ion chain is given by 

Ho = Y^ ^cr| + ^^o^AaLaa,. + ^), (1) 

j av 

where cr| is the Pauli spin matrix at the jth. ion site and 
the second term is the phonon Hamiltonian Hph with the 
phonon creation operator of the normal mode v along the 
three spatial directions a G X, Y, Z . The notation x, y, z 
refers to the pseudospin orientation in the Bloch sphere. 
The ath spatial component of the j/th ion displacement 
operator (5i?" is related to the ath phonon normal mode 
amplitude (unit norm eigenvector of the dynamical ma- 
trix) 6"'' and the ath phonon creation and annihilation 

operator via bk^ = Y.u \J 2mJz\-°-ay + aj^^] with M 
the mass of the ion and ujav the normal-mode frequency. 

A laser-ion interaction is imposed to create a spin- 
dependent force on the ions by using bichromatic laser 
beams to couple these clock states to a third state via 
stimulated Raman transitions Effectively, this pro- 
cess is equivalent to an off-resonant laser coupling to 
the two clock states by a small frequency detuning ji 
determined by the frequency difference of the bichro- 
matic lasers. The ions are crystallized along the easy 
axis (Z-axis) of the trap with hard axes in the X- and 



y-directions where the transverse phonons lie. Then cou- 
pling the Raman lasers in the transverse direction min- 
imizes effects of ion heating and allows for an identi- 
cal spin axis for each ion [IF. By accurate control of 
the locked phases of the blue detuned and red detuned 
lasers with similar Rabi frequencies, an effective laser-ion 
Hamiltonian [5^, along the spin direction can be 
engineered in the Lamb-Dicke limit Afc|(5Rj(t)| <C 1 : 

N 

ULiit) = 17j"Ak • (5Rj(t)(7j sin(/it), (2) 

i=i 

in which the effective Rabi frequency is generated 
by one effective blue-detuned beam and one red-detuned 
beam simultaneously (refer to Appendix A for details). 

In experiments, one uses adiabatic quantum state evo- 
lution to evolve the ground state from an easily prepared 
state to the desired complex quantum state that will be 
studied. For spin models generated in an ion trap, it is 
easy to create a fully polarized ferromagnetic state along 
the z-direction via optical pumping, and then apply a 
spin rotation (for instance, with a pulsed laser) to re- 
orient the ferromagnetic state in any direction (usually 
chosen to be the y-direction). Then, if one introduces 
a Hamiltonian with a magnetic field in the direction of 
the polarized state, it is in the ground state of the sys- 
tem. By slowly reducing the magnitude of the field and 
turning on the spin Hamiltonian of interest, one can adi- 
abatically reach the ground state of the Hamiltonian (at 
least in principle). Hence, one has an additional Zeeman 
term with a spatially uniform time-dependent effective 
magnetic field B(f) coupled to the different Pauli spin 
matrices as 

N 

nB{t)^J2^{t)-a, (3) 

(the magnetic field has units of energy here, since we 
absorbed factors of the effective magnetic moment into 
the definition of B, and it can also be expressed in 
units of frequency if we use units with h — 1). Note 
that we are using an unconventional sign for the cou- 
pling to the magnetic field, since this is the sign con- 
vention often used in the ion-trapping community. In 
this case, the ground state of the magnetic field Hamil- 
tonian has the spins aligned opposite to that of the field, 
while the highest-energy state has them aligned with 
the field. The magnetic field B(i) is made in the y- 
direction by directly driving a resonant radio-frequency 
field with frequency ujq between the two hyperfine states 
to implement the spin flips |23] or by indirect Raman 
coupling through a third state to effectively couple the 
two hyperfine states [52]. The full Hamiltonian is then 

n{t) ^nph + nLiit) + nB{t). 
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B. Factorization of the time evolution operator 

We solve for the quantum dynamics of this time- 
dependent Hamiltonian by calculating the evolu- 
tion operator as a time-ordered product U(t,to) = 
7texp[— i J^^dt"H{t')/?i] and operating it on the initial 
quantum state \ip{to)). For the adiabatic evolution of 
the ground state, we start our system in a state with 
the spins aligned along the magnetic field and the sys- 
tem cooled down so that there are no phonons at time 
to- IV'(^o)) = I tyty ■ ■ ■ ty) \0)ph- This means we will 
be following the highest excited spin state of the system, 
as described in more detail below. While it is also possi- 
ble to examine incoherent effects due to thermal phonons 
present at the start of the simulation, we do not do that 
here, and instead focus solely on intrinsic phonon cre- 
ation due to the applied spin-dependent force. 

In time-dependent perturbation theory, one rewrites 
the evolution operator in the interaction picture with re- 
spect to the time-independent part of the Hamiltonian. 
This procedure produces an exact factorized evolution 
operator 

i7(i,to) = 6-^^"^ (7Ki,io) (4) 

which is the first step in our factorization procedure 
[the first factor is called the phonon evolution oper- 
ator Uph{t,to).] [Note that Hph is time independent 
and it is multiplied by the factor {t — Iq) in the expo- 
nent.] The second factor is the evolution operator in the 
interaction picture, which satisfies an equation of mo- 
tion given by ihdUi{t,to)/dt = [Vi{t) + HB{t)]Ui{t,to), 
with Vi{t) = exp[zHp,i(i - <o)/?i]^L7(i)exp[-i'Hp/t(t - 
to)/h] since \HB{t)^^~iph] — 0. The only difference be- 
tween HLi{t) and V7(t) is that the phonon operators 
are replaced by their interaction picture values: a^i, — >■ 
aai, exp[-iuJa„{t - to)] and a^^ exp[iwQ,^(i - to)]- 

We now work to factorize the evolution operator fur- 
ther. Motivated by the classic problem on driven har- 
monic oscillators [24] , we factorize the interaction picture 
evolution operator via Ui{t, to) = exp[—iWi{t)/h]U{t, to) 
with Wi{t) defined by Wi(t) = ]l^dt'Vi{t') (we call 
the factor on the left the phonon-spin evolution opera- 
tor Uph-sp = exp[— zVF/(t)/?i] and the one on the right 
is the remaining evolution operator). The key step in 
this derivation is that the multiple commutator satisfies 
[[Wi{t),Vi(t%Vi(t")] = 0. This fact greatly simplifies 
the analysis below. 

The equation of motion for the remaining evolution 
operator U{t,to) satisfies 

in-u{t,to)^n{t)u{t,to). (5) 

in which the operator 'H(i) is given by the expression 
n{t) = e^^'^'^-ifidt +'HB{t) + Vi{t)]e-^'^'^'\ (6) 



The operator 'H(t) can then be expanded order by order 

as 

e^^^(*Vz(i)e-^^^(*) - Vi{t)+'^[Wi{t),Vi{t)] (7) 

N 

et^^W^sWe-^^^**) - ^{B(t).aj-f (8) 

'-[Wi{t)Mt)-o,] + ...] 

" ' 

Residual terms , 

e*'^^W^;^9te-sl^^(*) = Vi{t) 

+ \'^[Wi{t),vm (9) 

where we used the facts that dWi{t)/dt = Vi{t) and 
[[Wi{t),Vi{t')],Vi{t")\ = 0. Explicit calculations then 
yield 

N N 
j — l v—l a 

X sin(/xt)CTj, (10) 

in which r]ai, = Sax^k" y^^h/2~Muj^ is the Lamb-Dicke 
parameter for the phonon mode av with the ath com- 
ponent of the laser momentum Sk". When the terms in 
Eq. ([s]) vanish, virtually excited phonons will be shown to 
play no role on the spin-state probabilities as a function 
of time, but in the presence of a transverse field, due to 
the noncommuting nature of quantum operators, phonon 
creation can significantly affect the spin-state probabili- 
ties. This fact has not been considered in detail before, 
and involves one of the most important results of our 
work, as detailed below. 



C. Ising spin models and M0lmer-S0rensen gate 
for vanishing transverse magnetic field 

With a vanishing transverse magnetic field, the Hamil- 
tonian H can be greatly reduced to the spin-only Hamil- 
tonian Hspin{t) = i^\Wi{t),Vi{t)\. Because the spin op- 
erators s| in Wi and Vj commute, one can exactly derive 
the following Ising spin Hamiltonian [TU [T3] 

N 

HspUt)^ E Jn'(f)^l^r- (11) 

Then the expression for the spin exchange interaction 
Jn' (t) is 



4 



ujxu COS 2^t — 2fi sin uJxi^t sin /it] , 



(12) 



which can be uniformly antiferromagnetic {Jjj' (t) > 0) or 
ferromagnetic {Jjj'{t) < 0) for the instantaneous ground 
state of the Hamiltonian Hgpin (t) when the laser detun- 
ing /i is detuned close to center of mass mode frequency 
LOcM- However, the interaction Jjj'(t) > can also be 
inhomogcneous and frustrated when the laser is detuned 
in between phonon modes with details depending on the 
properties of the nearby phonon modes ujxi>, ^f^^ 

The M0lmer-S0rensen gate [T7] was originally proposed 
to disentangle phonon effects from the spins in ion-trap 
quantum computing. It was discovered that because the 
phonons arc harmonic, one could operate on the spins 
in such a way that the phonon state is unmodified after 
the gate operation (irrespective of the initial population 
of phonons). But one needs to keep in mind that this 
gate has no transverse field present, which can modify 
it because the transverse field operator does not com- 
mute with the Ising Hamiltonian. We begin our discus- 
sion by assuming that the laser is closely detuned to the 
transverse center of mass mode with angular frequency 
^CM (Ia'~'^ca/| ^ A*) and the addressing laser intensity 
for each ion is uniform and moderate (J7j = f2 <C /x). 
In this situation, the time dependent term with the fre- 
quency /J, -I- LOcM in the interaction Vj can be neglected. 
Therefore, the interaction Vi{t) and the operator Wi{t), 
which are proportional to the collective spin operator 



sj, can be reduced to the following forms 



Vi{t) = ^^^54e-(— ^)*acA/ + /i.c.],(13) 



VK/(i) 



2VN{ncM - Ai) 
X [(1 - e-*("^^'^-^)*)acM + h.c' 



(14) 



where rycM is the Lamb-Dicke parameter for the center 
of mass mode and to = is chosen. 

There are two parameter regimes where phonon ef- 
fects disappear. In the weak-coupling regime rjcM^ ^ 
/i — LUcM, the operator Wi{t) almost vanishes and the 
time evolution operator Uj is solely determined by the 
spin-only Hamiltonian because the phonon dynamics are 
adiabatically eliminated. Any extra phonon state re- 
distribution takes a long time to be experimentally ob- 
servable and therefore phonon effects are under control. 
Outside the weak-coupling regime, one can also prevent 
phonon effects by preparing spin states determined by the 
spin Hamiltonian Hspm (t) at a particular waiting time in- 
terval t = T. The idea is to choose the time interval T 
such that the operator is periodic with integer K cycles 



so that Wi{KT) = Wi{0). Therefore, the initial phonon 
state at the start of the simulation will be revived at the 
time intervals {ujcm ~ fJ-)T = K x 277 as can be clearly 
seen from Eq. ( |14[ ) when the we start with the phonon 
ground state |0) for example, but it is generally true for 
any occupancy of the phonon states. 



D. Effective spin Hamiltonian with a transverse 
field 



In the presence of a transverse magnetic field, the ideas 
of the M0lmer-S0rensen gate are modified. While it is 
tempting to claim that the residual spin-phonon terms in 
the magnetic field are irrelevant in the Lamb-Dicke limit 
Vai' ^ Ij it is difficult to quantify this if the residual 
terms are relevant in the presence of the time dependent 
magnetic field B(t) • aj which can be large in magnitude. 
In fact, phonon effects often modify the time evolution of 
the spin states when a transverse field is present. How- 
ever, one can say that in cases where the integral of the 
field over time is small (which occurs when the field is 
small, or when it is rapidly ramped to zero) or when the 
B field lies along the x-direction only (or vanishes), the 
residual terms in Eq. ^ are irrelevant. For large de- 
tuning (weak-coupling regime), where |W7(t)|/?i ^ 1 or 
rjxu^j <C wxi/ — the residual terms are always higher- 
order perturbations with respect to the leading trans- 
verse magnetic field term B(t) • aj in the course of the 
quantum simulation. On those occasions, one can also 
consider the system as described by only the quantum 
Ising spin model in a transverse magnetic field. In gen- 
eral though, we need to determine how large the residual 
terms are, which often can only be done with numerical 
calculations. We illustrate this for a number of different 
cases below. The residual terms are 

-Hresit) = exp[iWi{t)/h]nB{t)exp[-iWi{t)/h] - HB{t). 

(15) 

The equation of motion for U can be written as 



in^UitM) = [^[Wl{t),Vl{t)]+HB{t)+Hres{t) 



X U{t,to)- 



(16) 



We perform the final factorization by writing U{t,tQ) = 
Uspinit,to)Uentit,to) where the spin-evolution opera- 
tor satisfies ihdtUspmit,tQ) = {i[Wj{t),Vi{t)]/2h + 
'HB{t)}Uspin{i,ta) and the entangled evolution operator 
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satisfies The spin evolution operator Uspin{t,to) becomes 

t,to). 

(17) 



d 



Uspin(t,to) = 7^ exp 



N 



N 



= Tt exp 



dt H spin (t ) 



(18) 



which is the third factor for the evolution operator of 
the Ising model in a transverse field and we define 
Hspin{t) in the exponent. The spin exchange terms 



Jjj'{t) = Jjj, + AJjj'{t) as given in Eq. (12 1 include a 
time-independent exchange interaction between two ions 
Jjj, and a time-dependent exchange interaction AJjji{t). 
The time-independent term can be thought of as the ef- 
fective static spin-spin Hamiltonian that is being simu- 
lated, while the time-dependent terms can be thought of 
as diabatic corrections, which are often small in current 
experimental set-ups, but need not be neglected. For 
simplicity, we set the initial time Iq = 0. 

The entanglement evolution operator Uent is a com- 
plicated object in general, but it simplifies when one 
can approximate the operator Uent as Uent ~ 1 for the 
special situations discussed at the end of the last sub- 
section. In general, this evolution operator involves a 
coupling of spins to phonons in all directions and has a 
very complicated time dependence. If one evaluates the 
first few terms of the series for the time-ordered product, 
one finds it involves multispin interactions, spin-phonon 
coupling, and spin-exchange interactions in all spatial di- 
rections. But the net weight of all of the terms is gov- 
erned by the integral of the magnetic field over time, so if 
that integral is small, then this factor will also be small. 
Therefore, the adiabatic elimination of phonons based 
on M0lmer-S0rensen gate [T7] can be justified only in the 
case of a vanishing transverse magnetic field. With a 
constant magnetic field, the entanglement between spins 
and phonons can be periodic so that phonon effects can 
continue to be nulled at integer multiples of the appropri- 
ate period [HJ |T7j . But such a procedure would be more 
complicated than the standard gate, and is not relevant 
for adiabatic state creation simulations, so we won't dis- 
cuss it further here. From a mathematical standpoint, 
because the entanglement evolution operator is on the 
far right of the factorization, it's main effect is to modify 
the state from an initial spin state in direct product with 
the phonon vacuum to a state that will typically involve 



some degree of entanglement between phonon and spin 
degrees of freedom. 

We can use this factorization to show that in cases 
where the spin-entanglement evolution operator can be 
approximated by the unit operator, then phonons have 
no observable effects on the probability of product states 
(regardless of the number of coherent phonons created 
during the simulation), so this result is similar in spirit 
to the original M0lmer-S0rensen gate, but is different 
because it holds in the presence of a transverse field 
and requires no special times for periodic variations to 
recur. To do this, we need one final identity. We further 
factorize the entangled phonon-spin evolution operator 
eyip[—iWi{t)/h] into the product exp[— i r3f^(i)aj,^] 
exp[-i rxAt)a:o.] exp[-(l/2) TxAt)T*xM 
with the spin operator defined to be r^^(f) = 
J2j 'yxji't)'^j ^'^d its complex conjugate is Txu{t), while 
the function 7xj(0 satisfies 7xj(^) = ^jVXvb'^j/ifJ-'^ — 
ujxu) ^ [cTcpl— iu!x^t}{iuixu sin + ^ cos fit) — /i]. 

At this stage, we have factorized the evolution opera- 
tor into four main terms, each term being an evolution 
operator evolving the system from time to to time t. We 
have explicit values for the first three factors, but the last 
term (the entanglement evolution operator) can be quite 
complicated; we have also described situations where the 
exponent of that term is small and can be neglected. In 
this case, the probabilities to observe any of the 2^ prod- 
uct states with a quantization axis along the Ising axis 

= I tx or Ij.) (g) I ta; or 1^) (g) for the N ionic 

spins) is unaffected by the presence of an arbitrary num- 
ber of real excited phonons (which are excited by the 
phonon-spin evolution operator). Using the fundamental 
axiom of quantum mechanics, the probability Pp(t) to 
observe a product spin state starting initially from 
the phonon ground state |0) and not measuring any of 
the final phonon states involves the trace over all possi- 
ble final phonon configurations 



p,sit) = E - E - E 



nxi=0 nxv 



{nxi..nx...nxN\U{t,tom ® |*(0))| = 



(19) 
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where |$(0)) is any initial spin state (it need not be a 
product state) and nxu denotes the number of phonons 
excited in the i/th mode in the X-direction. The opera- 
tor in the matrix element entangles the phonons and the 
spins, so we evaluate the matrix element in two steps: 
(1) first we evaluate the phonon part of the operator ex- 
pectation value, and then (2) we evaluate the spin part. 
Note that since the pure phonon factor of the phonon 
evolution operator exp[—iHpht/h] is a phase factor, it 
has no effect on the probabilities when evaluated in the 
phonon number operator basis, so we can drop that fac- 
tor. Next, the term exp[—i^^Txi^it)ax^'] gives 1 when 
operating on the phonon vacuum to the right, so it can be 
dropped. We are thus left with three factors in the evo- 
lution operator. One involves exponentials of the phonon 
creation operator multiplied by spin operators (and is es- 
sentially a coherent-state excitation for the phonons with 
the average phonon excitation number determined by the 
spin state being measured) , one involves products of spin 
operators that resulted from the factorization of the cou- 
pled phonon-spin evolution operator factor, and one is 
the pure spin evolution factor Uspin- The two remaining 
factors that appear on the left involve only crj spin oper- 
ators, and hence the product state basis is an eigenbasis 
for those operators. This fact allows us to directly evalu- 
ate the expression in Eq. (19). We expand the evolution 



of the initial state at time t in terms of the product-state 

basis \^{t)) = Usprnit,to - 0)|$(0)) = E^' C/3' (t) |^') 

with 1/3') denoting each of the 2^ product state basis 
vectors and c^'(i) is a (complex) number. Using the fact 
that the product states satisfy the eigenvalue equation 
crjl/?') = m^p,\f3') with eigenvalues m^^, = -1-1 (for | t^.)) 
or —1 (for I -Ix)), we arrive at the expression for the prob- 
ability 



Pp{t) ^ \c0{t)\' exp 



1^3 j' 



N 



E n 



1 



jj' 



nxi---nxN=0 1^=1 

We used the matrix element 

{nxi-nxi.:nxN\e^^ r^^(*)'»x. |0) 

\nxi+---nxv + --nxN 



(20) 



i-^y 



^/nxil-nxJ.-nxN^- 



XN 



(21) 

in the derivation. The summations become exponentials, 
which exactly cancel the remaining exponential term and 
finally yield Pp{t) = \cp{t)\'^ , which is what we would 
have found if we evaluated the evolution of the spins 
using just the spin evolution operator Uspin and ignor- 
ing the phonons altogether. Hence, the coherently ex- 
cited phonons have no observable effects on the proba- 
bility of product states for the transverse-field quantum 



Ising model when we can neglect the entanglement evo- 
lution operator. If we do not measure the probability of 
product states, then the terms from the coupled spin- 
phonon evolution operator remain spin operators, and 
one can show that the probabilities are changed by the 
phonons. In other words, it is because the spin-phonon 
evolution operator is diagonal in the product space basis 
for phonons and spins that allows us to disentangle the 
phonon and spin dynamics. In cases where this cannot 
be done, we expect the phonon and spin dynamics to re- 
main entangled. In other words, phonons have observable 
effects on any spin measurements which introduces spin 
operators away from the Ising quantization axis such as 
most entanglement witness operator measurements. Fi- 
nally, we may ask what does the entanglement evolution 
operator do to the system? It is difficult to find any 
simple analytic estimates of the effect of this term, but it 
acts on the initial state which has the spins aligned along 
or opposite to the magnetic field and has no phonons. 
During the evolution of that operator, new terms will be 
created which involve entanglement of spin states with 
states that have created phonons. If the amplitude of 
those extra terms is small, they will not have a large ef- 
fect, but if it is not, then one has no other recourse but to 
examine the full problem numerically, which is what we 
do next. First we examine a perturbation-theory treat- 
ment, and then we consider the full numerical evolution 
of the system. 



E. Diabatic effects from time-dependent spin-spin 
interaction 



One may have noticed that the spin evolution opera- 
tor was not the evolution of a static Ising spin model. 
There were additional time-dependent factors in the evo- 
lution operator which arose from the additional time de- 
pendence of the exchange operators that was inherited 
by the phonons when they were "adiabatically" removed 
from the problem. In this section, we use adiabatic per- 
.turbation theory (reviewed in Appendix B) to analyze 
the effect of those extra time-dependent terms on the spin 
evolution of the system. In an adiabatic quantum simu- 
lation, one initially prepares the system in a certain pure 
state \n{tif)) of the initial Hamiltonian H{to) with the 
occupation a„(to) = 1 and the probability amplitudes in 
all other states vanishing [am{to) = 0]. Thereafter, the 
^ probability amplitudes to be excited into the other states 
can be approximated by 



dt 



E.,n{t') - E„{t') 



for later times, as long as the transition amplitudes 
|cKm(i)| are much smaller than one during the time evo- 
lution. This is the main expression we will use to eval- 
uate the diabatic effects due to the time-dependent ex- 
change interactions Jjj/{t). Here Em,n{t') are the instan- 
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taneous eigenstates of the spin Hamiltonian Hspin{t) with 
a static exchange interaction Jj,j'{t) = Jjj, and 9m,n(t') 
are the corresponding dynamic phases given by the inte- 



grals Or 



and we assume there are 



no degeneracies in the instantaneous spectrum as a func- 
tion of time. 

Let us briefly describe the experimental protocol for a 
typical trapped ion quantum simulator restricted to the 
spin-only Hamiltonian Hspin{t) = His(t) + HB{t) defined 
in Eq. (18). The system is initially prepared in a spin- 
polarized state I ty •■• ty ■•■ ty) along (or opposite to) the 
direction of the transverse magnetic field B(t) = B{t)y 
by optical pumping followed by a 7r/2 spin rotation. The 
spin-only Hamiltonian is then turned on with a maximum 
effective transverse magnetic field |B| — ^ \Jjj'{t)\ 
followed by an exponential ramping down of the mag- 
netic field to a final value Bme~*^^ at time t (r is the 
exponential ramping time constant for the decay of the 
magnetic field). After evolving to time t, the projection 
of the spin states along the a;-axis of the Ising Hamilto- 
nian is taken to find the probability to be in a particular 
spin state at time t (in actual experiments another tt/2 
pulse is applied to rotate the x-axis to the 2-axis where 
the measurement is made). 

If the system is perfectly adiabatic during the evolu- 
tion, the outcome of the quantum state would be the 
highest excited state of the Ising Hamiltonian Hjs{t) if 
the simulation starts out in the highest excited state 
of the magnetic field Hamiltonian Bm time 
t = io = 0, which corresponds to the spins aligned 
along the y-axis. This procedure is theoretically iden- 
tical to the ground state passage of the spin polarized 
state I ... 4,y ... ]^y) to the ground state of the negative 
of the Ising Hamiltonian —Hjs{t) with the system Hamil- 
tonian being modified as H{t) —H{t) [TTj. In a typical 
trapped ion quantum simulator, the frequency /i is suf- 
ficiently far from any phonon frequencies such that the 
condition rjxv^i <SC \^ — uxv\ holds to avoid the heat- 
ing of the system away from the initial phonon vacuum 
state during the simulation. In addition, the maximum 
magnetic field strength is much larger than the time in- 
dependent exchange interactions \Bm\ ^ \ to ensure 
the system initially starts in an eigenstate of the initial 
Hamiltonian. To optimize the adiabaticity of the simula- 
tion, the ramping time constant t for the magnetic field 
has to be chosen to be much greater than the largest 
characteristic time scale of the system, which is shown 
below to be the minimum of the inverse of the frequen- 
cies \n-Ulx„\,Bra/h. 

We now discuss the effects of the time-dependent ex- 
change interaction AJjj'{t). For concreteness, we will 
follow the highest energy state, starting from the spin 
state aligned along the direction of the magnetic field. 
Starting with the expression for the transition probabil- 
ity amplitude am (t) in Eq- ( 22 ) , we find the dominant 
diabatic transition is to the state with the minimum en- 
ergy difference AE — E„i{to) — £'„(io) with the initial 
spin polarized state | ty ••■ ty ••• ty)) assuming the matrix 



element in the numerator does not depend too strongly 
on (to|, which is true when Hsit) ^ His{t). At the ini- 
tial time t — 0, where the Ising couplings Jjj'{t — 0) can 
be shown to always vanish, all of the spin states with one 
spin fiipped along the y-axis are degenerate. This de- 
generacy will be broken by the Ising Hamiltonian Hjs at 
finite time t > t^. Due to spin-spin interaction in Hjs{t), 
the states along the y axis of the Bloch sphere called 
|m), have nonzero matrix components {m\dt' Hjs{t')\n) 
with the lowest energy gap AE — 2i3(to) with respect 
to the initial spin state \n) — \ ... ... ^y) . 

To approximately evaluate the transition amplitude 
am from the initial state to the two spin-flipped states 
|m), we do not actually need to know the state |m). The 
only relevant information we need is that it is one of the 
two spin-flipped states which tells us what the denomi- 
nator is. Hence, we can approximate 



{m{t')\dt'His{t')\n{t')) {m\dt'Hj,{t')\n) 



Emit') ~ E^it') 



2B{t) 



Using similar reasoning, we approximate A9{t') as 



Ae{t') 



dt- 



. 2B„ 



-2ujbt{1 



(23) 



(24) 



in which ujb = is the magnetic angular frequency. 
The operator dfHisit') = j' 9t'Jj,j'{t')'^j'^j' consists 
of modes with frequencies ujxv + <^Xiy ~ ^-^d 2/i with 
the time derivative dfJjj'it') given by 



N 



Xv 



sin 2/it' — cos wxzyi'sin ^t' 



^Xi 



smujxvt' cos /ii' 



-N 2 _,,2 — —i^^^Xuil-- — )suy{uxu - ti)t . 



2 ^ 



U> 



Xv 



^Xv 



(25) 

The last approximate expression is derived by keeping 
the contribution from the slow mode wxv — M ^-nd drop- 
ping the high frequency modes loxv + Ai, and 2/z because 
the detuning > is closely detuned to certain phonon 
modes in the quantum simulation. 

As a consequence, the probability amplitude Q;^^(t) 
induced by a single phonon mode is given by the expres- 
sion 



gEi 

n' 



m\a^al\n)m. 



(26) 

The function f{t) = 2i J^^ dt' sin {ujx^ - f^)t'e"^'>^*">+*-^ 
can be approximated in experiments [when ujxii — /i is 
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much larger than ljb at slow ramping {luxi^ — f^)T ^ 1] 
as 



fit) 



i{^Xu - 



(27) 



We therefore reach the conclusion that the probability 
amplitude a^^l ^, (t) is given by 



1 ^ ^f),,17,vry2 fof "^6^^ 

is I ,\ \ / \ X X \ \^ J -J 'yi-V J J 



3,3' 



X [1 - cos(a;x,y - fJ.)t] 



(28) 



We note that diabatic effects manifested in due 
to time-dependent Ising couplings grow exponentially in 
time as e*/"^ signifying that the theory is only accurate 
for short times. To suppress the diabatic effects, the cri- 
terion that has to hold for all phonon modes v is 



" 3,3' 



\n) 



< 1. 



(29) 

Based on this expression, when the laser is closely 
detuned to one of the phonon resonance frequencies 
oJXu^ the transition probability between states caused by 
Jj,j'{t') becomes large (diabatic). In addition, a stronger 
magnetic field is required to suppress the diabatic transi- 
tions with smaller detuning /i. This is supported by the 
following numerical discussion in section III C. Notice 
that the above expression should be a reasonable estima- 
tion as long as the condition B{t) ^ max| Jj j/ {t)\ = Jmax 
applies at time t after the beginning of the quantum sim- 
ulation. We can estimate the maximal time tc for which 
B{t) ^ max| Jjj'(i)| = Jmax holds. The cutoff tc is set 
by B{tc) = Jmax where Jmax is the absolute value for 
the maximum exchange interaction Jjji (t) between the 
spins. As a result, the cut-off time tc is proportional to 
the ramping time constant t with a logarithmic factor 
given approximately by 



tc = T In 



Jn 



(30) 



In the parameter regime where Bm/ Jmax 3> 1, in which 
our theory holds, tc can be extended somewhat beyond 
the ramping time constant r. Based on our numerical 
discussion, the diabatic effects are largest when the mag- 
netic field is ramped through the transition from para- 
magnetic state to other targeted spin states, which is 
also accompanied by larger phonon creations due to the 
shrinkage of the spin gap near the transition. 



III. NUMERICAL RESULTS 

In this section, we focus on showing the circumstances 
where quantum emulators can or cannot be described by 



the transverse-field Ising model with high fidelity. Since 
our goal is to understand under what circumstances the 
effect of the phonons is small, we consider different cases 
for the time-evolution of the system including various 
detunings and initial transverse magnetic field strengths. 

To isolate different effects, we compare two spin-only 
models in the presence of the ramping magnetic fields 
with the theoretically exact spin-phonon model based on 
numerical diagonization. The first is the ideal spin model 
which considers the evolution of the system with a static 
Ising model (spin-exchange coefhcients are the time- 
averaged exchange coefficients) and a time-dependent 
magnetic field. While one might think this is a purely 
adiabatic model, it has some diabatic effects, since the 
fully polarized state is not generically the ground state 
of the Ising plus magnetic field Hamiltonian because the 
(static) Ising exchange interactions are nonzero at the 
initial time. Hence, one can invoke a sudden approxi- 
mation to the system initially, and find that the initial 
state is a superposition of different energy eigenstates. 
In addition, the magnetic field varies in time and hence 
can cause additional diabatic effects due to its derivative 
with respect to time. 

The second is the effective spin model which involves, 
essentially, evolution of the spin system according to the 
spin evolution operator only in Eq ( |18[ ). Hence it has the 
static Ising Hamiltonian, the time-dependent Ising inter- 
actions and the time-varying transverse magnetic field. 
This model can have its Schrodinger equation solved in 
a spin-basis only, since all phonon effects are neglected 
except virtual phonon excitations. 

The third model is the exact spin-phonon model, where 
we evolve the system according to the original spin- 
phonon Hamiltonian expanded in the Lamb-Dicke limit 
[Eq. ([2])]. The only approximation used in this last model 
is the cutoff for the phonons. The strategy we use is to 
numerically integrate the Schrodinger equation using a 
direct product basis which involves a spin state in direct 
product with a phonon state. We do this because the 
Hamiltonian only connects states that differ by plus or 
minus one phonon number, and hence is block sparse in 
this basis. The spin states are chosen to include all pos- 
sible Ising spin states for the number of ions in the trap. 
The phonon basis is chosen to have a cutoff of a maximal 
phonon excitation. The maximal cutoff is always chosen 
to be larger than the average occupancy of the phonons in 
each normal mode of the ion chain. Of course, we expect 
more phonons to be excited into the phonon modes clos- 
est to the beatnote frequency of the lasers, so the cutoffs 
that are chosen will vary from one normal mode to an- 
other. For example, we often find we can set the phonon 
cutoff to be one for some of the phonon modes far from 
the driving frequency of the spin-dependent force. 

To facilitate our discussion, we define the root-mean 
square average Jrms of the fully connected Ising interac- 
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tion for N ions as 



N{N -1) 



(31) 



in which the static Ising interaction Jj^j' is given by the 
static term in Eq. (12 1 and the integer indices both 
range from 1 to TV. 



A. Symmetry and experimental protocols 

Let us discuss the symmetry of the spin-only system 
first, which is relevant for the exact diagonization of the 
spin-only Hamiltonian. There is one spatial inversion 
symmetry (Rj — ^ — R-:/) in the ion chain, since the equi- 
librium ion positions are distributed symmetrically about 
the origin in the trap and all phonon modes involve sym- 
metric or antisymmetric displacements of corresponding 
ion positions. There is also a spin reflection symmetry 
(cr^ — > — cr^, cr^ iT^, and cr^ — > "O"^) in the spin-only 
models with a transverse magnetic field B{t) J^j fj- This 
spin-reflection symmetry preserves all commutation rela- 
tions of the spin operators and leaves the Hamiltonian 
invariant. Therefore, there are four symmetry sectors 
for the eigenstates of the spin model (even space, even 
spin; even space, odd spin; odd space, even spin; and odd 
space, odd spin). 

If the static Ising couplings are all negative (positive) , 
the spin ground state is ferromagnetic (antiferromag- 
netic) and the highest spin eigenstate is the opposite, 
namely antiferromagnetic (ferromagnetic). We will focus 
on a detuning to the blue of the center-of-mass phonon 
mode. In this case, all spin exchange couplings are posi- 
tive and the ground state is antiferromagnetic, while the 
highest excited state is ferromagnetic. We will examine 
the adiabatic state evolution of the highest eigenstate. 
With all the respected discrete symmetries, we can con- 
struct the symmetric and antisymmetric ferromagnetic 
states of the spin-only Hamiltonian as ^(| fa; ••• tx 

)+\ix ■■■ ix)) 01- ^{\tx---'tx)-\ix--- ix)) which is 

in the (even, even) or (even, odd) sectors, respectively. 

The experimental protocol is to prepare the system 
initially in a spin polarized state | ' ' ' ty): [which is 
the highest eigenstate of the transverse magnetic field 
B{t)y], by optical pumping and a coherent spin rotation 
and then to gradually turn off the magnetic field with 
an exponential ramp B{t) = iJ^e"*/"^ while keeping the 
spin-dependent laser force in the cc-direction on during 
simulation time through stimulated Raman transitions 
between the spin states. According to adiabatic evolu- 
tion, if the quantum state is initially prepared in the high- 
est eigenstate of the field-only Hamiltonian B{t) J2j fjj 
the outcome of the quantum simulation will adiabati- 
cally follow the corresponding highest eigenstate of the 
Hamiltonian (Ising spin Hamiltonian), if there are no 
level crossings (which does not occur in this system). 



In the case with positive static Ising coupling, the fer- 
romagnetic highest energy eigenstate is the symmetrical 
ferromagnetic entangled state (the so-called GHZ state) 
j^{\tx---tx) + \ix--- ix)), when B>0. 

There are two intrinsic errors which can impede the 
quantum simulation in trapped ions. The first is diabatic 
effects which occur primarily when either parameters in 
the Hamiltonian are changed too rapidly in time, or when 
energy gaps in the instantaneous eigenvalue spectrum 
become to small. The second is the error induced by 
phonons in the presence of time-dependent transverse 
magnetic fields. For example, the phonon-spin Hamil- 
tonian does not have spin-reflection symmetry because it 
is linear in the operators, and hence the spin-phonon 
interaction breaks this Z2 symmetry. One consequence of 
this is to couple the symmetric and antisymmetric ferro- 
magnetic states which is likely to reduce the spin entan- 
glement of the GHZ state. (Other errors such as phonon 
decoherence effects due to spontaneous emission are not 
considered here.) 



B. Probability and GHZ state entanglement 
measurements 

Current experiments use atomic cycling transitions to 
measure the spin state of the ion (which clock state the 
ion is in), and do not measure the phonons excited in 
the system. Hence, the experimental observables are the 
probability Pp{t) of a spin-polarized state after tracing 
out phonons in the tensor product of the spin-phonon 
Hilbert space \spin)<Si\phonon) , as mentioned above when 
discussing Eq. ( [T9| . If one performs rotations about the 
Bloch sphere prior to making the measurement of the 
probabilities, then one can also measure a number of dif- 
ferent spin-entanglement witness operators. 

A spin-entanglement witness operator (for a target en- 
tanglement state) [5S] is a mathematically constructed 
observable that has a negative expectation value when 
the system is entangled. No witness operator can mea- 
sure general entanglement, but instead a witness operator 
is constructed to measure a specific type of spin entan- 
glement. For example, the witness operator VFghz for an 
TV-ion chain can be constructed as [25l 



Wgbz - 3/ - 2 



O +11 9 



k=2 



(32) 



with the stabilizing spin operators expressed in terms of 
the Pauli spin operators by 



N 

k=l 



'k-l<^k- 



(33) 



Notice that the target spin polarized state in this paper 
is along the Ising x axis in the Bloch sphere instead of 
the z axis. Therefore, we modified the original expres- 
sion [2 5) to our problem by the transformation cr^ — 
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and (T^ — >■ a^. Based on the above construction, GHZ 
state entanglement measurements can be detected by the 
observable (MKghz) = Tr{pWG}iz) with the density ma- 
trix p constructed by pure states or mixed states during 
the quantum simulation. For a perfect GHZ state entan- 
glement, one can show that the entanglement witness op- 
erator satisfies (Wghz) = ~1 (refer to Appendix C). Any 
deviation from perfect GHZ entanglement would lead to 
a value greater than —1. Note that this is one of the few 
cases of a witness operator where the degree of entangle- 
ment is correlated with the magnitude of the expectation 
value of the witness operator. 



C. Transition to the ferromagnetic state when 
detuned blue of the center-of-mass phonon mode 

The systems we consider range from iV = 4 to = 8 
which is far from the thermodynamic limit. The quantum 
phase transition (QPT) due to the discontinuity of the 
ground-state wave function in the thermodynamic limit 
N ^ oo only manifests itself as a state avoiding crossing 
in the energy spectrum, which is adiabatically connected 
to the QPT at large A^. The system parameters and the 
higher set of transverse phonon modes, which belongs to 
the higher branch of two transverse motional degrees of 
freedom, for different numbers of ions A^ are summarized 
in Table I. The trapping parameters are given by the as- 
pect ratio and the CM mode frequency ujcai along the 
transverse (tight) axis. The axial (easy) trapping fre- 
quency is given by the product of the aspect ratio and 
the CM mode frequency ujcm- The choice of these pa- 
rameters comes from trap parameters and typical oper- 
ating regimes of the ion-trap experiment at the Univer- 
sity of Maryland. Most results are robust with moderate 
changes of parameters and our choices do not intimate 
that fine tuning of parameters is needed to achieve the 
results we show. 

TABLE I: Parameter set I 



Aspect ratio 


0.2092 


ujcm/^tv 


4.63975 MHZ 


Rabi frequency ^l/2n 


369.7 KHZ 


Lamb Dicke parameter tjcm 


0.06 


No. of ions N 


4 


6 


8 




1.0 


1.0 


1.0 




0.9788 


0.9788 


0.9788 


Transverse 


0.9482 


0.9481 


0.9480 


phonon mode 


0.9088 


0.9083 


0.9079 


frequencies 




0.8589 


0.8581 


(in units of ujcm) 




0.7988 


0.7974 
0.7236 
0.6324 
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FIG. 1: (Color online.) Quantum simulation for different 
models (blue line: ideal spin model, green line: effective 
spin model, red line: exact spin-phonon model). The aver- 
age phonon occupation number n{t) is plotted with a solid 
black line in the right panels (e-h). The right panels also 
show the expectation value of the GHZ entanglement witness 
operator. The detuning p is always chosen to be 1.0095tJCM, 
just slightly blue of the CM phonon mode, while the expo- 
nential ramping time constant for the magnetic field is fixed 
at r = 8 X 10~^ msec. The initial transverse magnetic field 
varies in the different panels as follows; (a) Bm = 23.64Jrms, 
(b) B„ = 47.28 J™s, (c) Bra = 94.56 J™,, and (d) Bm = 
189. 12 Jrms. The average static Ising coupling (in units of an- 
gular frequency) is given by Jrms ~ 1.5017 x 10~*a;c]\f (for 
the chosen ion trap parameters in Table I). 



In Fig. [T] the left panel shows the probability to be all 



spins up, or all spins down [PFM(t) = \P\ 



|P^^^^...^^(t)p] and the right panel shows the correspond- 
ing entanglement witness operator expectation value 
{Wghz} for the GHZ state with four ions (A = 4), 
positive (blue) detuning (/z > locm) and positive Ising 
coupling (J > 0). In this case, the highest excited state 
is ferromagnetic and the center of mass (CM) phonon 
mode produces nearly uniform Ising couplings between 
each spin (more uniform the closer the detuning is to the 
CM frequency). 

We begin by examining the effect of the magnitude 
of the initial magnetic field on the adiabatic state evo- 
lution. In theory, we want the magnetic field as large 
as possible, because the larger it is, the better the ini- 
tial state is an eigenstate of the system. But this limit 
must be balanced against two competing issues: first, we 
(t)P -|- must complete the simulation in a finite time period and 
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FIG. 2: (Color online.) Detuning dependence /i of the simu- 
lation (the initial magnetic field is Bm = 3.6 x W'^locm and 
the exponential ramping time constant t — 8 x W~^msec 
are both fixed) for fi — I.OOSSujcm , Jrms = 3.765 x 
10-*wcA/ (a), = 1.0095a;cM, Jrms = 1.5017 x IQ-^wcm (b), 
/i = 1.019a;cAf, Jrms = 7.4734 x W'^locm (c), and fi = 
im8ujcM,Jrms = 3.7019 X W'^iJcM (d) . The other pan- 
els (e-h) are the corresponding GHZ witness operator and 
the average phonon number. 



we do not want parameters in the Hamiltonian to change 
too rapidly (to prevent large diabatic effects) and second, 
the experimental setup can only create a magnetic field 
of some finite maximum value. Hence, we must always 
live with a finite initial magnetic field. If the field is too 
low, then the sudden approximation tells us we will have 
larger admixtures of excited states in the system, which 
will give rise to oscillations in the measured probabilities, 
and will reduce the final probability to end up in the fer- 
romagnetic state. This is because the initial state is not 
close to the initial eigenstate of the system. If the field 
magnitude is much less than the spin coupling, then the 
system never evolves to the ground state, but remains in 
the state ordered along the magnetic field. 

The set of curves in panel (a) of Fig. [T] is at an inter- 
mediate value of the magnetic field. Let us focus first 
on the ideal spin model. One can see that the proba- 
bility has initial oscillations, which continue as the sys- 
tem evolves, but it produces a nearly pure ferromagnetic 
state probability by the end of the simulation. The cor- 
responding witness operator shows that we achieve excel- 
lent entanglement as well. If we add the time-dependent 



exchange interactions as described by the effective spin 
models we see the probabilities generate much larger di- 
abatic oscillations, and the final probability is strongly 
suppressed. In addition, the witness operator gives an 
expectation value shifted significantly up from —1. Fi- 
nally, when we consider the exact coupled spin-phonon 
model, we find the oscillations lie in between those of 
the ideal and effective spin models, but the probability 
for the ferromagnetic state is nearly one. The witness 
operator, on the other hand, shows reduced entangle- 
ment, similar to the effective spin model. This is easier 
to understand for the spin-phonon model, because the 
presence of phonons breaks the Z2 symmetry of the spin 
reflection, and hence the classification of states as even 
or odd is lost, so those states are mixed together by the 
Hamiltonian reducing the entanglement. Furthermore, 
comparing to the phonon occupancy, we see deviations 
between the red and green curve start when the phonon 
occupation is still quite low, and the oscillations in the 
witness operator are strongly correlated with oscillations 
in the phonon number (because the phonon and spin are 
coupled in a correlated fashion, of course). As the field 
is increased for fixed ramping time (panels b through d) , 
we see that the amplitude of the oscillations is reduced 
for the ideal spin model and for the exact spin-phonon 
model, but the witness operator for the exact solution 
always shows reduced entanglement. This behavior pro- 
vides one of our main conclusions about the effect of the 
phonons — namely, they act to provide a decoherence of 
the spin system, which reduces oscillations and improves 
the adiabatic nature of the evolution of the probabilities, 
but reduces the overall entanglement of the spin state. 

In Fig. [2] we examine the detuning dependence of the 
simulator in the case of a moderate magnetic field. We 
start with the detuning very close to the CM mode, and 
then move further and further away in panels (a) through 
(d). As expected, when we are far detuned from the 
phonon line (panel d), no phonons are created, and all of 
the three models agree essentially perfectly for the proba- 
bility and the entanglement. As we move in closer to the 
phonon line, we start to generate phonons and we also 
generate diabatic oscillations in the probability data, but 
surprisingly, we continue to see a very close correspon- 
dence of the exact spin-phonon model to the ideal spin 
model, while the effective spin model describes the sys- 
tem in an increasingly poorer fashion. The entanglement 
gets sharply reduced, as well, so much so that for the 
smallest detuning, the entanglement witness operator is 
effectively zero for the exact spin-phonon model. As be- 
fore the ideal spin model is incorrect in predicting the 
entanglement, and the effective spin model is much bet- 
ter at approximating the entanglement. 

In Fig. [3j we also examine the detuning dependence, 
but now in the case of a large magnetic field. Here the 

exact spin-phonon model and the ideal spin model give 
virtually identical results for the ferromagnetic probabil- 
ity and the effective spin model fails once the detuning 
gets too close to the phonon line. The entanglement. 
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FIG. 3: (Color online.) Detuning dependence of the sim- 
ulation with a strong initial magnetic field (Bm = 2.88 x 
W~^UJCM and r = 8x 10~^ msec) for /i — 1.0057locm , Jrms = 
2.5076 X IQ-^ujcM (a), n = 1.0095a;cM, J™^ = 1.5017 x 
10-*wcM (b), ^l = 1.0190a;cAf, J™. = 7.4734 x IQ-^cjcm (c), 
and A* = 1.0380a;cM, Jrms = 3.7019 x IQ-'^tJCAf (d). Panels 
(e-h) are the corresponding GHZ witness operator and the 
average phonon number. 



however, continues to behave as before, with the ideal 
spin model always predicting large entanglement and the 
exact spin-phonon model showing much reduced entan- 
glement. Note that the trend of reducing the ferromag- 
netic probability and the entanglement as the detuning 
moves away from the phonon line is arising from the fact 
that the exchange parameters are getting smaller, so the 
field is evolving the system too rapidly to be in the adi- 
abatic limit and we have diabatic effects which do not 
allow us to reach the ferromagnetic state during the time 
of the simulation. 

If we go to the case of near resonant driving of the 
system, as in Fig. |4j we see that we create significant 
numbers of phonons, and we eventually see that the exact 
spin-phonon model does deviate from the ideal spin model 
in the probabilities. As expected, we also see a very sharp 
reduction of the entanglement of the exact evolution as 
compared to the ideal spin model. 

So far, all of our calculations were for a small chain 
with = 4. Now we look at the size dependence of the 
ion chain by examining the behavior as we increase the 
number N from 5 to 8 in Fig. [5] This case is with a mod- 
erate initial magnetic field, but driving fairly close to res- 
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FIG. 4: (Color online.) Quantum simulation with detuning 
near resonance with the CM mode and with a large initial 
magnetic field {Bm = 2.88 x 10~^ljcm and r 
for 11 = 1.0038CJCM, Jrms = 3.765 x lO"*- 
1.00095a;cM, Jrms = 1.5 X lO'^iJCM (b). Panels (c-d) are the 
corresponding GHZ witness operator and the average phonon 
number. 



onance with the CM mode. In order to properly compare 
the different systems, we adjust the detuning to yield an 
average spin exchange that is approximately equal for the 
different number of ions. These results clearly show the 
main themes we have been exploring. The exact spin- 
phonon model and the ideal spin model exhibit nearly 
identical evolution of the probability, with the deviations 
starting to become clear only for the largest chain. Nev- 
ertheless, the entanglement of the system is sharply re- 
duced for all systems, with the reduction occurring at 
about the same time in the simulation for all cases. In 
terms of CM mode phonon excitations, we observe the 
trend of increasing phonon excitations from far below one 
phonon toward more than one phonon as the system size 
increases. The strong correlation between phonon exci- 
tations and the deviation between the exact spin-phonon 
model and the ideal spin model is still clearly manifested. 



D. Evolution to the kink phase in the highest 
excited state 



Let us look at the highest excited state on ion chains 
with even numbers of ions and a detuning that lies 
between the highest (CM) mode and the second (tilt) 
mode. We examine this scenario here to see the effects 
of phonons on this transition. 

In Fig. |6j we show predictions for experimental traces 
when the detuning is between these two phonon lines. 
The red detuning is chosen such that it is possible to ex- 
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FIG. 5: (Color online.) Number dependence of quantum sim- 
ulation for intermediate initial magnetic field and detuning 
chosen according to the text. In the left panel subplots, we 
plot the ferromagnetic probability as a function of time. In 
the right panel subplots, we plot the time dependence of the 
GHZ state witness operator (VKghz). Red solid lines are re- 
sults for the exact spin-phonon model and blue solid lines are 
for ideal spin model. The black lines indicate the correspond- 
ing coherent CM phonon occupation driven by the lasers. 
The initial magnetic field strength Bm = 3.6 x W~'^locm 
and exponential ramping time constant t — 8 x 10^^ msec 
are chosen with varying parameters given by the following 
cases: (a) N ^5,fi = 1.0076a;cM, Jrms = 1.5013 x lO-"*- 



(h) N = = 1.0063uJCM,Jrms = 1.5011 X 10'Wm, (c) 
N = 7,fi= im54ujcM,Jrms = 1.5006 X W^ujcM, and (d) 
N ^ 8,fi = 1.0048t^cM, Jrms = 1.5008 x IO^Wm- Panels 
(e-h) are the corresponding GHZ witness operator and the 
average phonon number. 



FIG. 6: (Color online.) Quantum simulation for the kink 
state in the presence of phonon excitations. The left panel 
shows the probability of the one-kink state Pk as a function 
of time for the ideal spin model (blue), effective spin model 
(green) and the exact spin-phonon model (red) for different 
cases. The right panels show the average phonon occupation 
number n^{t) for the CM mode (black lines) with angular 
frequency lucm = 2tt x 4.6398 MHZ and the tilt mode (blue 
lines) with angular frequency ojt = 0.9788cjcm. We fix the 
exponential ramping time constant t = 8 x 10~^ msec for dif- 
ferent magnetic field and detunning as: (a) /i = 0.9905tJcM, 
= 3.6 X IQ-^iJCM, (b) fi = 0.9810a;cM, Bm = 3.6 x 
IQ-^OJCM, (c) ^ = 0.9905a;cM, Bm = 2.88 x 10"^cjca/, and 
(d) n = 0.9810a;cM, -Bm = 2.88 x W^ujcm- Panels (e-h) 
are the corresponding GHZ witness operator and the average 
phonon number. (System parameters are chosen according to 
Table I.) 



cite two transverse phonon modes simultaneously by lo- 
cating the laser detuning between the CM mode on an ion 
chain with iV = 4 ions ujcm ~ 27r x 4.64 MHz and the tilt 
mode with frequency lut — 0.9788ujcai ■ The spin state 
simulation ideally should evolve from the highest excited 
state I ty ■■• ty ■■■ ty) along the transverse magnetic 
field axis to the symmetric kink state (or antiphase state) 
^(1 tajt^^J-a^ix) + I ixixtxtx)) along the Ising axis if the 
system is perfectly adiabatic. In those cases, long range 
spin-spin couplings are negative and dominant over adja- 
cent spin couplings; therefore, the highest excited state is 
the symmetric kink state (with ferromagnetic couplings, 
the highest excited state will be antiferromagnetic) . By 



detuning from the CM mode tdcM to the tilt mode ujt, 
the adjacent spin-spin couplings change sign from neg- 
ative to positive and long range spin-spin couplings are 
still negative with enlarged magnitude. We have checked 
that this eigenstate is generic for even number of ions 
(up to at least = 8) with the highest eigenstate always 
being the antiphase (kink) state. 

In cases (a) and (b), with the same magnetic field 
ramping, we observe a larger probability for the kink 
phase Pk in panel (b) (further detuned from the CM 
mode) for the ideal spin model. Similarly, we also ob- 
serve a larger Pk (higher fidelity) in case (d) than in 
case (c). Lower fidelity in general is correlated to the 
fact that the spin excitation gap is smaller (with respect 
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to the exponential ramping time constant r). This im- 
plies that the frustration of the Ising couplings favors a 
larger spin gap for intermode detuning than that for the 
close detuning case to the blue of the CM mode. By exact 
numerical diagonization for ideal spin models, we found 
it is indeed the case and the second highest eigenstate 
has very different character. In cases (a) and (c), the 
second highest eigenstate is the symmetric antiferromag- 
netic state ^(| txixtxix) + \ ixtxixtx}) at the end of 
the simulation with a spin gap « 27r x 28.359KHZ 
. However, for cases (b) and (d), the second highest 
eigenstate is the kink state ^{\ txlxlxlx) + \ ia^txta^t^: 
) + I irEi2;ixt2;) + I ti;txt2;ii)) with a larger spin gap 
As w 27r X 107.40KHZ. 

As far as phonon effects are concerned, we notice that 
the large discrepancies between the ideal spin model and 
the exact spin-phonon calculations for Pk are strongly 
correlated to the large number of phonon excitations 
(flu > 1) for either phonon mode. By comparing case 
(a) with case (c) [or comparing case (b) with case (d)], 
we notice that large phonon excitation often occurs dur- 
ing the part of the evolution when the magnetic field 
strength is larger. Unlike the blue detuned situation in 
subsection C, a small detuning from the CM mode does 
not show high fidelity of the final spin state (kink order 
in this case) at long times. Instead, the probability for 
observing the kink state is significantly degraded by the 
presence of the phonons. Hence, the creation of phonons 
does not always help the system appear as if it is simu- 
lating a static spin Hamiltonian. 

By greatly reducing the Rabi frequency to suppress the 
phonon excitations, we observe good agreements between 
ideal spin models and exact spin phonon models in case 
(a) as well as the enhancement of the spin probability 
Pk- For larger ion number case {N = 6), we find the 
diabatic effects completely ruin the simulation by tuning 
accessible system parameters in Table I to an order of 
magnitude larger or smaller. Even if one can suppress 
the phonon effects in this case, we will still typically have 
large diabatic effects due to the small energy gaps. 



E. Kink transition in the ground state 

Following the ground state evolution with an odd num- 
ber of ions and a laser detuning between the second and 
third highest transverse phonon modes, there is a sharp 
phase transition |26j separating ferromagnetic and one- 
kink spin ground states for the ideal spin models. There- 
fore, it is quite relevant to examine the feasibility for 
experimentally observing this transition by going beyond 
an adiabatic evolution of the system and treating phonon 
effects and diabatic effects exactly. The crossover be- 
comes sharp very quickly with the number of ions due to 
the fully-connected nature of the spin models. A typical 
adiabatic phase diagram for the spin only Hamiltonian 
with iV = 3 and = 5 is shown in Fig. [7] The order 
parameter Pp — Pk for the transition is defined as the 



difference between the probabilities of the two degenerate 
ferromagnetic states Pp and the four degenerate one-kink 
states Pk ■ The one-kink states for the A = 3 case are 
given by the four degenerate states | \xixix), \ ix\x\x), 
I \rx\rx\x), and I \x\x\rx)- For A = 5, the four one- 
kink states I ^-Axixixix) ,1 ixix^x^x^x), I ixVxix^x^x), 

and I tajtxta^iKii) are also degenerate. It is clear that 
the crossover between the two phases is sharper (smaller 
spin gap) for A^ = 5 than = 3 as shown by the much 
sharper contrast in color. In addition, the ferromagnetic 
phase area moves closer to the leftmost phonon mode and 
shrinks in size as the number of ions increases. As the 
magnetic field gets smaller, the transition between the 
ferromagnetic and kink phases becomes very sharp as N 
increases [ST . The numerical discussion in this section is 
based on the system parameters in Table II. 



TABLE II: Parameter set II 



Aspect ratio 


0.1 




4.63975 MHZ 


Rabi frequency n/27r 


3.697 KHZ 


Lamb Dicke parameter ricM 


0.06 


No. of ions N 


3 5 


Transverse 


1.0 1.0 


phonon mode frequencies 


0.9950 0.9950 


(in units of to cm) 


0.9879 0.9879 




0.9789 




0.9683 



In Fig. [Sj we numerically map out the time depen- 
dence of the probability Pp — Pk for ideal spin mod- 
els. The nexponentially ramped magnetic field B{t) is 
chosen with different initial values Bm (scaled by Jrms 
as determined by the detuning /x and the Rabi angular 
frequency fi). The value of the Rabi angular frequency 

= O.OluJcM is chosen so that it is safely within the weak 
field regime "q^fli, <Si {\ujxi' — /^l) near the central region 
of the phase diagram for all phonon modes. The total 
simulation time T is chosen so that it is proportional to 
the inverse of Jrms/'^TT. We select the exponential ramp- 
ing time constant r for the exponential reduction of the 
magnetic field B{t) to be one-fifth of the experimental 
simulation time (r ~ 0.2T). By comparing (a) to (b) [or 
(c) to (d)], we observe that diabatic effects are greatly 
suppressed when the exponential ramping time constant 
T is large enough so that the transition to the closest 
excited state in energy is negligibly small. This effect 
shows up as much deeper colors dictating the order pa- 
rameter Pp — Pk in the ferromagnetic states and the 
kink states when B{t) approaches zero on the vertical 
axis, as illustrated in subplots (b) and (d). The diabatic 
effects also show up clearly as a slow oscillation in the 
probabilities Pp — Pk at larger B{t)/Jrms > 1 before 
the simulation ends along the vertical axis in subplots 
(a) and (c). We also notice some fast background os- 
cillations in Pp — Pk covering the entire phase diagram 



15 



#4 



N=5 




I 



□ .989 [].gg 0.991 D.992C.993 0.994 



D.989 0.991 0.992 0,9930.994 

^'/"coM 



FIG. 7: (Color online.) Adiabatic phase diagram of the trans- 
verse field Ising model for (a) = 3 (left panel) and (b) 
N = 5 (right panel). The horizontal axis is the laser detun- 
ing scaled by the transverse CM mode frequency ujcm- The 
vertical axis is the transverse magnetic field B scaled by the 
root-mean-square average of Ising coupling Jrms- The blue 
area represents the one-kink phase and the red area indicates 
the ferromagnetic phase. The range of the detuning ^ (in 
units of loom) is shown between the second phonon mode 
and the third phonon mode. The value of the order param- 
eter Pf — Pk ( varying from —1 to +1) is described by the 
color scale to the right of the figure. 
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FIG. 8: (Color online.) Ferromagnetic to kink phase diagram 
for N = 3 calculated for the ideal spin model with diabatic 
effects included. The horizontal axis is the laser detuning fj, 
scaled by transverse CM frequency ujcm = 2n x 4.6398MIIZ. 
The vertical axis is the instantaneous transverse magnetic 
field B{t) scaled by the root-mean-square average of the spin 
couplings Jrms (uotc the range changes for different panels). 
The Rabi angular frequency 57 and the dimensionless Lamb- 
Dicke parameter for the center of mass mode rjcM are se- 
lected to be = O.OlcJoAf and rjcM = 0.06 respectively, 
(a) r = 0.2r,r = 0.614 x ( J™,/27r)-\ = 5J™s, (b) 
T = 0.2T, T = 6.14 X {Jrms/2n)'\T = 0.2r, B„ = 5J™s, (c) 
r = 0.2r,T = 0.614 x (J™/27r)-\r = 0.2T, = lOJrms, 
and (d): r = 0.2T,T = 6.14 x (Jrms /2n)-\ B^ = lOJ™^. 



in cases (c) and (d). This effect is due to the fact that 
the time derivative of the dynamic phase 6m{t) — 6'„(t) 
between (to, n) states is roughly stationary in time (as 
analyzed with adiabatic perturbation theory). For short 
exponential ramping time constants r, one cannot see 
the noticeable interference pattern between these states 
because of a random phase cancelation along the path of 
the state evolution in time. In panels (b) and (d) the 
main difference is a reduction of the period of the back- 
ground interference pattern, which is shorter in panel (b) 
(larger magnetic field). 

What are the phonon effects on the corresponding fer- 
romagnetic to kink phase diagram? We show our calcu- 
lations for the = 3 case in Fig. [9j Phonon creation 
has serious effects on the phase diagram. In the cases 
with fast ramping time constants [(a) and (c)], the fer- 
romagnetic states are destabilized and appear only with 
small probability. For slow ramping time constants [cases 
(b) and (d)], the FM domain disappears near the leftmost 
phonon mode due to large phonon creation as the phonon 
is being more resonantly driven. But the kink state do- 
main reduces only slightly near the rightmost phonon 
mode. As a consequence, phonons restrict the available 
parameter space to observe the FM to kink phase dia- 
gram but do not rule out the possibility of observing the 
phase as long as the exponential ramping time constant 



T is long enough. In the current numerical simulation we 
show, the exponetial ramping time constant r is roughly 
on the order of a few milliseconds (close to feasibility in 
current experiments). One may suspect that phonons 
can ruin the stability of the FM state when the number 
of ions scales up because the FM domain shrinks in size 
and moves closer to the leftmost phonon mode, and if we 
are too close to the phonon mode, phonon creation ruins 
the chance to see the FM state. One can try to increase 
the experimental simulation time and reduce the Rabi 
frequency, but doing this too much eventually runs into 
coherence issues or problems from spontaneous emission. 

In Fig. [To] we show the N = 5 case for the ideal spin 
model. The behavior of the FM-kink phase diagram is 
similar to the A^ = 3 case in Fig. |8] except the bound- 
ary of FM states and kink states is shifted toward lower 
detunings fi. As a consequence, the FM state domain 
(deep red area) occupies the region where the detuning 
is close to the leftmost phonon mode and the kink state 
domain (deep blue area) occupies most of the detuning 
region. However, the phase diagram for low magnetic 
field B{t) — 7> is very close to the adiabatic phase dia- 
gram (see the right panel of Fig. [?]) when the initial mag- 
netic field Bm is large and the ramping time constant t 
is long as shown in Fig. 10 (b) and Fig. 10 (d)] except for 
the background interference patterns that were described 
above. One also notices that there is much less diabatic 
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FIG. 9: (Color online.) Ferromagnetic to kink phase dia- 
gram for N — 3 with the same cases as in Fig[8] The phase 
diagrams are calculated with the exact spin-phonon Hamil- 
tonian. The horizontal axis is the laser detuning ^ scaled 
by transverse CM mode frequency ujcm- The vertical axis 
is the instantaneous transverse magnetic field B(t) scaled by 
the root-mean-square average of spin-spin coupling Jrms- (a): 
T = Q.2T,T = 0.614 X ( J™,/27r)-\ B„ = 5J™,, (b): r = 
0.2r,r = 6.14 X (J™,/27r)-\r = 0.2T,B^ = 5J™,, (c): 
r = 0.2r,r = 0.614 X {Jrms/2n)'\T = 0.2T,B^ = lOJ™,, 
and (d): r = 0.2T,T = 6.14 x ( J™^/27r)-\ B„ 
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FIG. 10: (Color online.) Ferromagnetic-kink phase dia- 
gram for A'^ = 5 as calculated for the ideal spin model. 
The horizontal axis is the laser detuning /i scaled by trans- 
verse mode trapping frequency ujcm — 2tt x 4.6398MIIZ. 
The vertical axis is the instantaneous transverse magnetic 
field B{t) scaled by the root- mean-square average of spin- 
spin coupling Jrms- The Rabi angular frequency Q and 
the dimensionless Lamb-Dicke parameter for the CM mode 
rjcM are SI = 0.0086a;cM and rjcM = 0.06, respectively, 
(a) r = 0.2r,r = 0.614 x ( J™,/27r)-\ B„ = 5J™s, (b) 

T = 0.2T, T = 6.14 X {Jrms/2TTy\T = 0.2r, Bm = 5 Jrms, (c) 
T = 0.2r,T = 0.614 X {Jrms/2TT)-\T = 0.2T, Bm = lOJrms, 

and (d) r = 0.2r,r = 6.14 x ( J™«/27r)-\ = lOJ™,. 



effects at low B{t) due to the fact that the smallest spin 
excitation gap is larger near the central area of the phase 
diagram. 

When we add phonon effects, we might expect the 
phase diagram to only deviate when we are detuned close 
to a phonon line, but the situation is much worse for 
= 5, as shown in Fig. 11 The kink phase (deep 
blue zone) exists for a wide range of ramping and on- 
set magnetic fields Bm as shown in all cases. However, 
the ferromagnetic domain (red zone) disappears even for 
slow ramping [like the exponential ramping time constant 
T w 30 msec in cases (c) and (d)]. This does not rule out 
the possibility of observing the FM phase for even longer 
ramping time constants t (or smaller Rabi angular fre- 
quency ri) but a ramping time constant t « 30 msec is 
already well beyond what is used in current ion-trap ex- 
periments where r is usually less than one millisecond. 
This problem gets worse for larger N, and already for 
N = 7 the FM-kink phase diagram appears to be impos- 
sible to observe. This arises in part due to the fact that 
the spin gap closes exponentially fast with the system 
size [15]. As a result, one needs to dramatically reduce 
the diabatic effects to see the transition. In addition, 
phonon effects also make it hard to see the transition 
by not allowing the detuning to move too close to cither 
phonon line, and thereby misses significant regions where 
the FM phase is stable. 
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FIG. 11: (Color online.) Exact FM-kink phase diagram for 
the spin-phonon Hamiltonian with N = 5 and the corre- 
sponding parameter set as in Fig. |10[ The horizontal axis 
is the laser detuning jj, scaled by the transverse CM fre- 
quency UJCM- The vertical axis is the instantaneous trans- 
verse magnetic field B{t) scaled by the root-mean-square 
average of spin-spin coupling Jrms- (a) r = 0.2r, T = 

0.614 X {.Jrms/2n)-\Bm = 5Jrms, (b) T = 0.2T, T = 

6.14X (J™/27r)-\r = 0-2T,Bm = 5J™,, (c) r = Q-2T,T = 
0.614 X {Jrms/2-Ky^,T = 0-2T, Bm = lOJ™,, and (d) r = 

0.2r,r = 6.14 X (Jrms/2ny\Bm = WJrms- 
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IV. DISCUSSION AND CONCLUSIONS 

We have examined a number of different issues related 
to the importance of phonons in analog quantum simu- 
lation of the transverse field Ising model. We show when 
the spin-phonon entanglement operator Uen can be ap- 
proximated by 1 (a longitudinal magnetic field along the 
Ising axis, or a vanishing transverse magnetic field), one 
can show that the phonons do not affect the probability 
to measure the spins in product states in the direction of 
the Ising interaction, but they can reduce the entangle- 
ment of the spin eigenstates. Surprisingly, in cases when 
the operator Uen cannot be approximated by 1, the ef- 
fect of the phonons is often to make the system look more 
like a static Ising spin Hamiltonian plus a time-varying 
transverse magnetic field. This result holds primarily 
when laser detuning is blue of the CM mode, and hence 
corresponds to a ferromagnetic case when one looks at 
the highest excited state. We emphasize that the com- 
mon belief based on the geometric phase gate, in which 
phonon effects can be suppressed by choosing the period 
to be the inverse of close detuning from a phonon mode 
due to periodic spin-phonon entanglement dynamics, is 
no longer valid in a finite decaying transverse magnetic 
field. 

Our work shows that one must consider phonon effects 
in most ion-trap spin simulator experiments especially 
when the spin-spin interaction is highly frustrated. In 
cases when laser is detuned blue of the transverse cen- 
ter of mass mode, phonons are beneficial to make the 
system look more and more like the static spin Hamil- 
tonian being emulated (at the expense of reducing spin- 
spin entanglement). In cases when spin interactions are 
frustrated with multiple phonon modes stimulated, the 
phonons can work to suppress the true adiabatic spin 
phases from having a high fidelity or even invalidate the 
spin phases completely. Generically the phonon effects 
beyond adiabatic elimination are remarkable when the 
detuning lies close to at least one of the phonon frequen- 
cies, and hence on average more than one phonon per 
mode is excited. 

In conclusion, large laser detuning is essential to sup- 
press phonon coherent population while it also causes the 
shrinkage of spin excitation gaps in the adiabatic spin 
simulation. Alternative adiabatic quantum simulation 
schemes which do not create noticeable phonon occupa- 
tion, while maintaining large spin excitation gaps would 
be desirable. 
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Appendix A: Spin-dependent force and effective 
transverse magnetic field 

We describe how different laser-ion interactions are em- 
ployed to ultimately simulate effective spin models. The 
spin-dependent dipole force along any direction of the 
equatorial plane of the Bloch sphere and an effective 
transverse magnetic field are created by using multiple 
laser beams and the optical-dipole interaction between 
the ions and phase-locked lasers. We start with a refer- 
ence Raman beam that has frequency and then su- 
perpose another perpendicular Raman beam (that has 
frequencies in a frequency comb lul + ujq + fi,ujL + ^^o ^ 
fj-jiOL +wo). The lasers use off-resonant Raman coupling 
through dipole-allowed excited states of the ion to gen- 
erate an effective spin-phonon interaction. Take the Yt- 
terbium ion ^''^Yb''" as an example: the qubit states are 
the clock states \F — 0,mi? — 0) and |F = l,mi? = 0), 
formed from the hyperfine states of the S valence electron 
and the spin one-half nucleus. These states have no linear 
Zeeman effect, and hence are less prone to background 
magnetic fields. 

The hyperfine state \F — 1, =0) is denoted as the 
spin-up state and the state \F = 0,mF = 0) is the spin- 
down state in the z-axis of the pseudospin Bloch sphere. 
The energy-level spacing is in the gigahertz range, so one 
could, in principle, directly make transitions that flip the 
spins from up to down and vice- versa by stimulated emis- 
sion and absorption processes acting on the hyperfine 
states. But, it is common to instead generate these spin- 
flip transitions via off-resonant Raman coupling to a third 
state to suppress incoherent spontaneous emission effects. 
To do this we need two laser beams with different fre- 
quencies which can be detuned away from the energy level 
spacing between the clock states and each of the frequen- 
cies are chosen to be far away from dipole allowed reso- 
nant transitions in each ion. We denote the two beams 
wavevectors and frequencies by (ki, wi), (k2, W2) respec- 
tively. By adiabatic elimination [55] of dipole-allowed 
excited states through the Raman procedure for ion j, 
one can write down the interaction for an ion as 



H 



LI 



i(Ak-Rj-Awt+<# 



H.C. 



(Al) 



where r2| is the effective Rabi frequency of the 
stimulated-Raman transition, hAk = ?i(ki — k2) and 
hAio = h{iOi — 0^2) > are the effective momentum and 
energy of the photons, respectively, (p is the controlled 
phase shift between the two laser beams, and the pseu- 
dospin flip operators are Sj' , S~ . 

The full Hamiltonian involves the sum of this term 
plus the clock state energy level difference hujo multiplied 
by the z-component spin operator. Now we go to the 
interaction picture with respect to the clock state energy 
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level difference Tuo, 



'0, 



HIj = exp[^a|t]i7i, exp[-^a|t] (A2) 

at time t. With the photon energy difference ?iAa; compa- 
rable to the clock-state energy splitting Tiwo, only terms 
with slow modes (rotating wave approximation) are kept 
and we arrive at the following Hamiltonian relevant for 
our discussion: 



ttRWA 



_J_q+ 

2 



,i(Ak-Rj(t)-<5 



H.C., (A3) 



in which the slow mode angular frequency is given by 
duj — Auj — ujQ and (j) is the static phase shift between 
the laser beams. The coupling of the reference Raman 
beam with photon frequency ujl and the blue-detuned 
photon with frequency wl -I- -t- /i (/x > 0) in the sec- 
ond beam leads to an effective blue-detuned beam with 
the frequency difference Aw — ujq + ii as given by the 
Hamiltonian Hi, 



N 

E 



■i(Ak-Rj(i 



^'^"^o'+'f"'^ +H.C.), (A4) 



where Hf, is the interaction with the blue detuned (Awh = 
wo + m) laser that has a beatnote frequency /i > and Ak 
is the wavevector difference of the two interfering laser 
beams that generate the Raman coupling, Rj (i) = R° -|- 

STljit) is the ion position operator with the equilibrium 
ion position at site j. Similarly, the coupling of the 
photon from the reference beam with the photon in the 
red-detuned beam with frequency ujl+ujq — fi leads to the 
effective red-detuned laser with the frequency difference 
Auj = luq — fi given by the Hamiltonian Hr 



N 

E 



+ gi(Ak-Rj(t 



)-Au.^t+4>r-) _^ H.C.). (A5) 



Employing a superposition of multiple frequency com- 
ponents and adiabatically eliminating the dipole allowed 
excited states allows one to show that the interaction 
of laser beams with ions consists of interactions between 
the reference beam and the other frequencies. As a 
result, after the summations in Eqs. (A4) and (A5), one 
arrives at the following expression: 



g-iAk-<5Rj(t) g-^+i<t>s 



(A6) 



in which hermiticity of the Rabi frequency is used, and the static phases are 05 = — Ak • R" — ((/if, -t- 4'r)/'2 and 

4>M = (0r — 06)/2. In the Lamb-Dicke limit, we have exp [lAk • STij] fa 1 + lAk • SRj, and the Hamiltonian HfJ^^ 
is reduced to 



H 



RWA 
LI 



nm cos(^t -f (j)M)S+e-''^^ 1 + i Ak • 5Rj {t) 



H.C. 



(A7) 



The first term only induces resonant carrier transitions 
in the pseudospin sector without coherent phonon exci- 
tations. The second term induces first (red or blue) side- 
band transitions with the change of one phonon occupa- 
tion number at each phonon mode as can be seen by re- 
placing the displacement operator 5'Rj (t) by phonon cre- 
ation and annihilation operators i ^i^av i ^ 

ajjj,). The spin-dependent force pointing along the az- 
imuthal angle (jjg in the equator of Bloch sphere is then 
derived from the phonon side bands as 

Hlf w nn'j cos(^i + (l)M)<jf''Ak ■ SKjit), (A8) 

in which the spin phase is given by 4>s = 4>s ~ |'^^ the 

relation S^e^^'^s + Sje^^'^s = is used, and the spin 

orientation is given by n = cos((/)g)x -I- sin(0g)y. 

The expression for the spin-dependent force can be jus- 
tified by keeping the phases 0^ = 0,(f>M = ""/S locked. 
Take a transverse phonon mode scheme for example, in 



which Ak • R° = 0, with Ak || (5Rj (t). The spin- 
orientation = can be locked along the x-axis in 
Bloch sphere when the phase difference <j>r — 7r/2,0;, = 
— 7r/2 is maintained. This can be achieved by passing 
the second beam through an acousto-optic modulator 
(AOM) maintaining the phase difference between the fre- 
quency components n^+LOQ, fi^+cjo+Z^, ^l+^o^ ^J■ to be 
out of phase. As one can tell from the dependence of the 
spin phase 4>s on Ak - Rq, it is not sensitive to transverse 
phonon excitations (coherent or thermal) in contrast to 
the sensitivity it has to the longitudinal phonon modes. 
This is why most state-of-art trapped ion quantum spin 
simulators couple to transverse phonon modes. 

One should note that there is a fast oscillating term 
in the transverse magnetic field —ft^'j siii{fit)a^ due to 
carrier transitions. This term causes very fast oscillations 
of low amplitude which are averaged over during the time 
of an experiment, so we neglect them here. 

Let us now consider how to generate a slow effective 
transverse magnetic field by using two continuous Raman 
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beams with frequencies a;^, and +ajo, with phase dif- 
ference 0, and wavevector difference Ak. Starting from 
Eq. (A6) but with a different effective Rabi frequency fl^ 
for the resonant beam with ii — 



Hfr = ^ 



^+^.Ak.[R^^R^(t)]+0_^^^^^^ (A9) 



we choose the lasers to be out of phase (0 = T7i'/2) so 
that the side-band terms vanish within the Lamb-Dicke 
expansion exp [iAk • S'R.j{t)] w 1 + i6k ■ S'Rj{t) . The 
effective transverse magnetic field can then be derived 
by direct substitution as 



] 3 ■ 



(AlO) 



in which the transverse magnetic field is given by 
Bj = ±Ml^ /2 when the phase shift cj) is equal to 
=F7r/2 and (t| = 2S^j is the Pauli spin operator (we 
will be working in a nontraditional Pauli spin matrix 
representation, where is diagonal, is real, and a' 
is imaginary). Hence, the transverse magnetic field By 
can have its amplitude changed as a function of time 
by adjusting the laser intensity in the mode that has its 
frequency equal to oj^ + with an AOM. 



Appendix B: Adiabatic perturbation theory 

The time-dependent Schrodinger equation for the evo- 
lution of the wave function \il>{t)) is (we drop the spin 
subscript on the Hamiltonian) 



(Bl) 



Since the Hamiltonian is always Hermitian, we intro- 
duce instantaneous eigenfunctions \n{t)) with the instan- 
taneous eigenenergies defined by 



H{t)\n{t))^E,,{t)\n{t)). 



(B2) 



The time-dependent wave function |'0(i)) can then be 
expanded in terms of the orthonormal eigenbasis \'n{t)) 
as 



m)) = Y.^n{t)\n(t)), 



after using the orthonormality relation {m{t)\n{t)) = 
5m,n for the instantaneous eigenfunctions. One can fur- 
ther relate the matrix elements {m{t)\dt\n{t)) to the ma- 
trix elements {m{t)\dtH {t)\n{t)) . Simply take the time 
derivative of Eq. (B2) and project onto the state |m(t)), 
to show 



{ni(t)\dMt)) 



{m{t)\dtH{t)\n{t)) 

En{t) — Em{t) 



(B5) 



for n ^ m [this derivation assumes the instantaneous 
energy spectrum has no states that are degenerate with 

En{t)]. 

In the adiabatic approximation, the transition matrix 
elements {m{t)\dtH {t)\n{t)) between different instanta- 
neous eigenstates are assumed to be so small they can 
be neglected. In this limit, the system simply follows 
the instantaneous eigenstates \n{t)) without transitions 
between different instantaneous eigenstates. In general, 
transitions between eigenstates |n(i)) should be consid- 
ered to determine the corrections to the adiabatic state 
evolution. Choose the gauge Cm{t) = am(t)e~*^'"*^*) with 
the phase 9m (t) — J^^ dt' ^ . The probability ampli- 
tude am{t) induced by the diabatic transitions can be 
solved by integration with respect to time in Eq. ( B4 ) as 



da,n{t) 
dt 



E 



{m{t)\dtH{t)\n{t)) ,[g^^t)-e„(t)] 



Em{t)-En{t) 



(B6) 



in which 9m{t) and 9n{t) are the dynamic phases accu- 
mulated by the states \m{t)) and \n{t)). 

In an adiabatic quantum simulation, one initially pre- 
pares the system in a certain pure state |n(io)) of the ini- 
tial Hamiltonian i?(io) with the occupation Q;„(to) = 1 
and the probability amplitudes in all other states vanish- 
ing [am{ta) = 0]. Therefore, the probability amplitudes 
to be excited into the other states can be approximated 

by 



{m{t')\dt.H{t')\n{t')) ,[e,„(0-9„(n] 
Emit') - E.^{t') ^ 



for later times, as long as the transition amplitudes 
|am(i)| are much smaller than one during the time evolu- 
tion. This is the main expression we will use to evaluate 
the diabatic effects due to the time-dependent exchange 



(B3) interactions Jjji(t). 



in which the coefficients C„(t) = {n{t)\tp{t)) are the time- 
dependent quantum amplitudes projected onto the in- 
stantaneous eigenbasis \n{t)). Therefore, the equation 
of motion for the expansion coefficients Cn(t) can be de- 
rived by direct substitution into the Schrodinger equation 
in Eq. (Bl) which becomes 



ih 



Cm{t) +Y,Cn{t){m{t)\dt\n{t))'j - Em{t)Cm{t), 

71 

(B4) 



Appendix C: GHZ state entanglement 

The observable for the measurement of GHZ state en- 
tanglement is given by 



{Wghz) = Tr {pWoHz} 



(CI) 



with the density matrix p constructed from either pure 
states or mixed states after tracing over the phonons. 
One can explicitly show that a pure GHZ state |GHZ) = 
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73 (I txtx ••• tx) + I ■■■ Ix)) leads to an entan- 

glement measure (M^ghz) to be —1. The entanglement 
measure (Wghz) is independent of the basis; therefore, it 
is natural to choose spin polarized states along the Ising 



X-axis in our discussion. Based on the GHZ state |GHZ), 
the corresponding density matrix pghz = | GHZ) (GHZ | 
in the Ising product state representation is given by 



Pghz 



~t~ I * * *)(4'a:4'X * * * I ~t~ I 4'x4'a: * * * 4'a:) (l^xl^x 

I 



(02) 



Therefore, the nonzero matrix elements in pcHZ arc with the stabilizing spin operators S'p^^'^ , S^^"^^ ex- 
among the fully ferromagnetic states and each has an pressed in terms of the Pauli spin operators by 
expectation value equal to 1/2. As a consequence, one 
only needs to know the Wqhz matrix elements among 
the FM states to evaluate the GHZ state entanglement 
(Wghz)- Following the definition of the witness operator 
Wghz 



N 



Sf^^^ = l[al S^''^-=at_,at, (C4) 



k=l 



Wghz = 31-2 



9 11 9 



the nonzero matrix elements for the global spin flipping 
(C3) operator ^^Zn amongst the N ion FM states is 



fe=2 



(tcctx ■ ■ ■ |<S'i'^^'~' I XxXx ■ ■ ■) — (4-x4-a; ' ' ' |S'p^^'* | tcctx ■•■) — !) 



(C5) 



in which the spin operator erf flips the spin state of i-th One can also derive the nonzero matrix elements for 
ion as tx) = | ix) and crfl ix>) = | tx) with I the the two-particle spin operator S^^^" as: 
unit operator for spin states. 



N 



N 



(txtx •••in ^r^" + ^1 ^x4x • • •) = iixix •••in ^k""^" + ^1 txtx • • •) = 2 



JV-1 



(06) 



fe=2 



fe=2 



As a conseqiience, the only nonzero matrix elements for 
the witness operator Wghz are the ones between two FM 
states which become 



following manipulations 



(txtx •••IWghzUxIx •••) 
= (txtx^^^|WGHz|txtx^^^) 
= -1. 



(07) 



Hence, the measurement (Wghz) for the pure GHZ state 
is characterized by the value —1 as can be seen by the 



(Wghz) = (tx • • • IpghzI tx • • •)(;x • • • |Wghz| tx • • •) 
+(tx • • • IpghzI tx • • •)(tx • • • IWghzI tx • • •) 
= i x(-l) + i X (-1) = -!. 

(08) 
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